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ABSTRACT 



During the annual Gray Whale migration from the 
Aleutians to Baja California, the mammals travel in coastal 
waters, thereby presenting an opportunity for the study of 
their sound spectral and source levels and vocabulary. How- 
ever, such measurements are distorted by surface and bottom 
reverberation. Using the theory of rough surface scattering, 
knowledge of the bottom impedance, and correlation techniques, 
it is possible to decompose the shallow water reverberation 
into the contributions from different paths . From this , the 
range, depth and the deverberated spectral source levels of 
the sounds of the mammal can be determined by use of only one 
hydrophone rather than the conventional three or four. The 
theory, deverberation programming, and experimental results 
are presented for a model of the whale's pulsed radiation in 
a laboratory model coastal environment. 
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I . INTRODUCTION 



Once yearly, the Gray Whale, Eschrichtius glaucus , 
migrates southward from the Aleutians and passes very close 
to the California coast in shallow water. During this 
migration, the sounds in the water close to the whales can 
easily be recorded; but, they may not be the true sounds 
produced by the whale. The whale normally produces an 
intermittent, short duration signal which in shallow water 
is received at the hydrophone via direct, surface reflected 
and bottom reflected paths. Since the path lengths are 
different, the signals arrive at the hydrophone at different 
times and they interfere. To obtain the true sounds produced 
by the whales, this interference, which is called reverbera- 
tion, must be eliminated. 

The purpose of this thesis is to model in the laboratory 
the whale sounds in the shallow water and to develop a 
technique to determine the range, eliminate the surface and 
bottom reverberation and calculate the spectral source levels 
as a function of time. 
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II. THEORY 



In a reverberant environment, the original signal can 
only be realized if the reverberation is removed. The 
method used to remove the reverberation, which is here 
called "deverberation , " assumes that the whale is a point 
source, the geometrical spreading is spherical, the water is 
isovelocity and the water depth is constant. 

Before the deverberation technique can be applied, the 
direct, surface reflected and bottom reflected path dis- 
tances must be known. Normally, this information is obtained 
by knowing the source position and calculating the path 
distances from the known geometry. Determining the hori- 
zontal source position is difficult and requires at least 
two directional or three omnidirectional hydrophones, accu- 
rately fixed with respect to each other at all times. To 
determine the depth requires at least one additional hydro- 
phone at a different depth. Generally, three or four 
hydrophones are deployed [Refs. 1 and 2]. With each addi- 
tional hydrophone, the complexity of the system is increased 
since each hydrophone requires its own amplifying and record- 
ing system. In shallow water, however, taking advantage of 
the surface and bottom reflections, one can use a single 
hydrophone and gather all the information necessary for the 
application of the deverberation technique. 

Consider the direct, surface scattered, and bottom 
scattered sounds received at only one hydrophone which 
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is deployed in shallow water as depicted in figure 1. It 
will be shown that when the differences between the arrival 
times for the three paths are known the three path distances 
can be calculated. Using the surface specularly scattered 
path, R , it is seen that 

R s - Re, 4 Rs 
R' s * =D’t5' 1 

where 

Tr»- Ch-0? 

(VUD^ L 

c" 2 ■ f R z - (.H - O ) 7 

5 Ch + d ) 1 L 

Substituting and rearranging terms gives 

(1) R s =TSTdT 

(2) R s-(trW k**™ 0 )* 

Using r , the time difference between the direct path 
arrival and the surface reflected path arrival and C, the 
mean speed of sound, gives 

c T s - R 5 + R" 5 - R 

and therefore 

x. 

(3) Ct s - + 4-UD)* - R 

Similarly using the bottom specularly scattered path 

Rb s Rb 4 Rb 

R'g - U-d ) 2 * B' 2 

r'4* = U - d) 2 + b“ 2 

o 
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FIGURE 1*. DEPLOYMENT OF ONE HYDROPHONE IN 

SHMLOW W&TER ENVIRONMENT 



where 



5 = U2.-H-D) 

r“- 



[ R 4 -Ch-D') 4 ]’* 

[r 4 - (.» -of ] 1 



Substituting and rearranging 

O') 



(4) = \z 2. - H ■ • D T I^ 2, "* ^ 2 + V\0 - - ?.D)j 

R » s U-H) r R » 4 4(? 4 +HD-lH-lD)l 

(5 > k B (2i -H -D) l K V J 



JL 

Z 



Using Tg, the time difference between the direct path 
arrival and the bottom reflected path arrival gives 



CT 






- R 



and, therefore, 

X. 

(6 ) Ct B + 4(l 2 +HD^H-^D)] - R 



Solving equations (3) and (6) simultaneously for R, the 
range of the source from the receiver and D, the depth of 
the source, produces 



(7) D = 

(8) R = 



( C 4 T?.'tp, -v 41 4, - - T c (c.r B ) 

4[HT B 4 T s U-h)] 

4nn- 

2CT S 



The equation for the range is left as a function of the water 
depth to facilitate its being programmed. Now that the range 
and depth are known, the other two path distances can be 
calculated. From equations (1) and (2) , the surface reflected 
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path distance is 



(9) R s = (R Z -V <*HDV 

and from equations (4) and (5) , the bottom reflected path 
distance is 

(1°) R b = \r 2 4 ^ WD - ZH - ID)] Z 

Therefore, when the time arrival differences for the differ- 
ent paths are known, equations (8) , (9) , and (10) can be 

used to determine the first three path distances. The paths 
for multiple reflections can be calculated directly from the 
known geometry, assuming specular scatter. 

For a transient signal, determination of the differential 
arrival times x^ and x^ can be realized by the use of an 
autocorrelation technique. The correlation function can be 
defined as [Ref, 3] 

(id r to -- e { O (.o - u] W u ) - u.1} 

where v is the time average, u is the mean, and E is the 
expected value of the received signal. This process is 
programmed using a digital summation 

n-X 

(12) R (O = S ( [ VV. w - all MUr W - a] 

with n representing the total number of samples in the record. 
Performing the autocorrelation on the reverberant signal at 
the one hydrophone gives peaks at delay times corresponding 
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to zero delay time, and the arrival delay times from the 
reflected signals. The peaks are realized only when the 
direct signal is delayed enough to correlate with the 
reflected signals. The computer program called AUTOPEAK 
performs the autocorrelation and then applies an envelope 
detection to determine the delay times for the peaks. These 
delay times are then used in equations (7) , (8) , (9) , and 
(10) to determine the range, depth, and reflected path dis- 
tances. Thereby, all the geometrical information necessary 
for the deverberation technique has been obtained. 

Computer programs have been developed for deverberation 
in either the frequency or time domain. 

The computer program designed to eliminate the reverbera- 
tions in the frequency domain is called DEVERB . For the time 
being, assume only one frequency component, whose amplitude 
and frequency are functions of time. For a signal with time- 
varying frequency and amplitude, we can describe the pressure 
at the receiver, due only to the direct path signal as [Ref. 
4] 

_ / S / V JwWi 

(i3) P D (O = c(i)e 

Then, taking into account spherical divergence, the spatial 
phase shift, and specular scattering from a Gaussian rough 
surface, the pressure at the hydrophone due to the surface 
scattered signal can be written as 

tu s 
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and the pressure due to the once bottom reflected signal is 



reflection coefficients, and and tt are the phase shifts 
due to the bottom and surface reflections, respectively. 

The surface reflection coefficient depends on the roughness 



of the surface. The exponent for specular scattering is 
[Ref. 5] 



where 0* is the R.M.S. wave height, A is the wavelength of 
the sound signal and 0g is the angle of incidence. Equations 
(14) and (15) are derived from previously received direct 
path pressures, corrected for path differences and phase 
shifts and represent the pressures due to the single reflected 
paths. The coherent sum of equations (13), (14) and (15) is 
the pressure sensed by the receiver when these three compo- 
nents are present, 

The signal processing is done digitally. Therefore, the 
continuous time dependence is replaced, whenever it occurs 
in the previous equations, by digital block numbers, indi- 
cated by the subscript index "k." Each block contains 
enough data samples of the incoming signal to give the desired 
spectral frequency resolution during a block duration which 
is small compared to the total duration of the time-varying 
signal. The frequency change of the signal within a block 




(R and e are the frequency- dependent pressure amplitude 



(16) 9 2 = cos e s 
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duration is assumed to be much smaller than the frequency 
resolution of the digital processing. The block duration is 



T 



_ NUMBER OF SAMPLES IN THE BLOCK 
SAMPLING FREQUENCY 



The time the block ends is related to the continuous time t 
by 

(17) t = TK; K = 1, 2, 3 

The index : 'i" is used for the spectral frequency compo- 
nent of the complex wave being analyzed. 

The equation for the source pressure at unit distance 
using the frequency deverberation correction is 



where l(K-N), l(K-M) and l(K-L) are unity factors with values 
l(K-N) =1 K > N 

= 0 otherwise 

l(K-M) =1 K > M 





0 otherwise 
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l(K-L) =1 K > L 

= 0 otherwise 



. th 



The pressure amplitude of the i frequency component in 
block K of the receiver reverberant signal is represented 
by ^ and its phase represented by <J>^ ^ . The deverberated 
pressure amplitude is represented by D„ . and its phase by 

IX , 1 

^ v . . The first term on the right hand side of the equation 
represents the received signal. The second term represents 
the correction due to a single specular scatter from the 
surface; the third represents a single bottom reflection 
correction, and the fourth represents the correction for a 
path which includes one surface and one bottom reflection. 

The equation can be expanded to include other multiple reflec- 
tions . 

The block indices are determined by 



(19) M = K - 

T 

(20) N = k - Zx 

T 

(21) L = K - -J SB 

T 

The output of the frequency deverberation program is a 
series of Spectra of the consecutive blocks and its Fourier 
transform is a time plot of the deverberated signal. 

The above procedure takes the signal from the time 
domain (the time series after A/D conversion) by Fourier 
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transform to the frequency domain, where the known frequency 
dependent reflection coefficients are easily applied, and 
then back to the time domain to verify the effectiveness of 
the process. 

When the reflection coefficients can be assumed to be 
frequency-independent, a simple point-by-point deverberation 
procedure can be applied in the time domain. The applicable 
temporal deverberation equation is 

< 2 2 ) D * = C . k *< 6 ETC - 

th 

represents the pressure amplitude for the K c sample. The 
other terms are similar to those in equation (18) . For low 
roughness surfaces g < 1 and the use of e over the 

appropriate frequencies will be a good approximation which 
is essentially independent of frequency. One advantage of 
the temporal deverberation technique is the relative freedom 
from restrictions of block size; the block size is determined 
only by the desired frequency resolution and the rate of 
change of frequency of the transient source. 
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III. PROCEDURE 



In order to model the Gray Whale's environment in 

Monterey Bay, a three meter cube "anechoic" water-filled 

tank was used. An artifical bottom made of hard rubber with 

£ 

an experimentally determined pc of approximately 2.4 x 10 

mks rayls was positioned one meter below the water surface. 

This type of bottom was chosen for its specific acoustic 

impedance since the bottom at the listening area in Monterey 

£ 

Bay has a pc of approximately 3 x 10 mks rayls. The depth 
of the bottom and relative placement of the source and 
receiver were determined in order to obtain realistic modeled 
delay times between the reflected signals and the direct 
signal. A 1.8 cm diameter spherical source was used because 
of its small size and its ability to transmit a signal with 
minimum distortion; but it also limited the minimum frequency 
to about 10 kHz. An FM up-sweep of varying widths and sweep 
rates was used to model one of the sounds produced by the 
Gray Whale. 

The equipment was connected as in figure 2- with the mas- 
ter unit pulse generator being used to synchronize the samp- 
ling frequency oscillator which determined the start and stop 
frequencies, sweep rate and pulse width. The pulse repeti- 
tion rate was determined by the master unit pulse generator. 
The signal was then amplified and sent to the source. The 
reverberent signal was received by an LC-10 hydrophone, and 
then amplified again and sent to the analog- to-digital 
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FIGURE 2= EQUIPMENT BLOCK DIAG RM\A 



converter. The A/D converter was triggered through the 
delay generator via the unit pulse generator and the sampling 
frequency oscillator. A sampling frequency of 320 kHz was 
used, The delay generator delays the trigger by the time 
required for the signal to be transmitted through the water 
and then this delayed pulse triggers the unit pulse generator. 
The unit pulse generator was used to determine the total on- 
time of the sampling frequency oscillator. The oscillator 
determined the samples per second taken by the converter 
during the oscillator's on-time. This complex equipment set- 
up was designed to allow the A/D converter to sample only 
during the time when the direct and reflected signals were 
present, thereby reducing the amount of computer time and 
storage required to process the data. Reflections from 
other surfaces in the tank were eliminated wherever possible 
in the sampling time by varying the pulse repetition rate 
of the master pulse generator or by varying the source and 
receiver placement. After the analog signal was changed by 
the converter to digital form, it was stored on cassette 
memory and then analyzed using the programs AUTOPEAK and 
DEVERB . 
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IV. DATA PROCESSING AND RESULTS 



The autocorrelation plots of reverberation for two 
different sweep widths are seen in figure 3 with the scale 
factor for the delay time equal to 3,125 y seconds. The 
plot of the 90 kHz sweep width shows a steeper slope of the 
envelope, thus a more clearly defined peak than the one with 
only a 10 kHz sweep width. This indicates, as expected, that 
as the difference between the upper and lower frequencies 
decreases, the correlation peaks become harder to determine. 

In the limit of only one frequency being present, there would 
be no peaks in the envelope. Figure 4 shows the range error 
(computer determination compared to direct measurement) versus 
the ratio of the upper frequency to the lower one. The graph 
indicates that for a ratio above 1.2 to 1 the frequency sweep 
of the signal is sufficient to get accurate time differences 
for the reflections and thereby to determine the source range 
and depth from the autocorrelation processes. 

The frequency deverberat ion program is designed to give 
a true spectrum of the signal by eliminating frequency depen- 
dent reverberations. Since the signal is time-varying, a 
small block of time is desirable to keep the change in the 
signal to a minimum during the block. The limiting factor 
to the minimum block size is the desired frequency resolu- 
tion. For a spectrum to accurately represent an instant of 
time rather than being a spectrum averaged over a length of 
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time, the frequency change during the block time should be 
much less than the frequency resolution. In addition, where 
possible, the block duration should be equal to, or a sub- 
multiple of |T - T^| in order to permit equation (18) to 
be applied. 

From figure 4, it is known that depth and range can be 
determined only when the frequency ratio of the chirp is 
greater than 1.2 to 1, Figure 5 shows the reverberant and 
deverberated signals for a chirp from 46 kHz to 54 kHz with 
a ratio of 1,17 to 1. The reverberant signal, top of figure 
5, was divided into blocks, equal to x fi - tg, and then trans- 
formed back into the time domain which is shown at the bottom 
of the figure, Some improvement can be seen but the expected 
slowly increasing frequency at the approximately constant 
amplitude has not been realized. It is believed that the 
inadequate deverberation is due to the fact that the frequency 
sweep during a block is approximately four- tenths of the 
frequency resolution. A slower sweep rate or a larger block 
duration would have cured this problem. However, a slower 
sweep rate for the model would have decreased the accuracy of 
the range determined by the autocorrelation; and a large block 
duration was precluded by the geometry which determined 
| Tg - Tg| , The temporal deverberation technique which is 
presented next did not suffer from these limitations . 

The result of applying the temporal deverberation program 
to a signal with a sweep width from 5 kHz to 95 kHz can be 
seen in figure 6 with the reverberant signal on the top and 
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FIGURE 5. 

FREQUENCY DEVERBERATION 
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FIGURE 6 



TEMPORAL DEVERBERATION 
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the deverberant signal below, Since the source was resonant 
at 64 kHz, the FM sweep grows in amplitude to 64 kHz and 
then decreases. The deverberated signal shows the frequency 
and amplitude modulation cleanly until the end of the direct 
signal. The reverberation after that time is due to scatter- 
ing from the side walls of the tank. Possible reverberation 
due to the fourth and later terms on the right hand side 
of equation (22) were excluded by limiting the duration of 
sampling by the computer. 



27 



V. CONCLUSIONS 



The work described in this thesis demonstrates the 
feasibility of obtaining a non-reverberant spectrum of a 
transient source in a reverberant environment. The technique 
includes calculation of the autocorrelation of the received 
signal to determine range and depth of the source and com- 
puter processing to correct for the surface and bottom reflec- 
tions , 

The autocorrelation function provides an accurate method 
for obtaining the range and depth of a source of transient 
sound in shallow water. The correlation technique can be 
performed for a chirp sound with ratio of upper to lower 
frequency of greater than 1.2 to 1. At least two reflections 
are required to obtain the depth and range of the source with 
respect to the receiver. 

Frequency and time deverberation programs which use the 
position data from the autocorrelation have been developed 
to eliminate the reverberations and, thereby, to obtain 
corrected spectra or corrected time plots. Either technique 
can be used; however, because the output of the computer is 
a time series, it is natural to apply temporal deverberation. 
This becomes very simple if the surface and bottom reflection 
coefficients are independent of frequency. 
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COMPUTER PROGRAM AUTOPEAK 



1 

2 

G 

10 

15 

20 

25 

30 

35 

40 

45 

50 

55 

100 

105 

107 

108 
110 
112 
112 
115 
120 
125 
127 
13C 
135 
140 
145 
150 
155 
16 C 
165 
1 7 C 
175 
180 
185 
190 
192 
205 
210 
215 
220 
225 
230 
235 
240 
245 
250 
255 
260 

265 
26 7 

266 
270 
275 
280 
285 
290 
295 
500 
503 
505 
5 1 C 

512 

513 

514 
51 5 
520 
525 

527 

528 
530 



V$ ( 4 8 1 . V 1 $ ( 6 ) , Z$( 1 ) 

A( 1000) ,6(1000) , Y( 1000) 
Z1S ( 1) 

P( 500, 



REM 
REN 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
REM 
DIM 
DIM 
DIM 
DIM 
REM 

Z4$ = ,, N" 
Z$="Y" 

REM 
REM 

GOSUB 500 
IF Z4$=" Y" 
REM 

GOSUB 1000 
REM 

GOSUB 1300 
REM 

GOSUB 1500 
REM 

GOSUB 2000 
IF Z3$="N" 
REM 

GOSUB 2500 
REM 
GOSUB 3000 
IF Z8=l GOTO 
REM 

4000 



DATE OF LAST CORRECTION: 1C/31/77 



* ALTOPEAK — 9/29/77 

* JEANIE SAVAGE, PROGRAMMER 

* SPECIFICATIONS BY RICK ECSTIAN 

* THIS PROGRAM PERFORMS AN AUTO- 

* CORRELATION ON SIGNAL DATA AND 

* THEN PICKS OUT THE TWO PEAK 

* VALUES. 

£ sjt ijs aj. # 3?$ a=s 



,Z2$(1),S(1),D(1),Z4$(1) 
1) , Z3$ ( 1 ) ,M( 50,2) ,R(1 ,2) 



***DRIVER ROUTINE*** 
PERFORM INITIALIZATION 



PROCEDURE 



'GOTO 140 
READ DATA 



FROM TAPE 



BUILD CTHER ARRAY 
PERFORM AUTO- CORRELATION 
DETERMINE INTERMEDIATE PEAK VALUES 



GOSUB 

REM 

GOSUB 

REM 

GOSUB 

REM 

GOSUB 

REM 

GOSUB 

PRI NT 

INPUT 

IF Z $= 

PRINT 



GOTO 185 

PRINT INTERMEDIATE 
DETERMINE TWO PEAK 



150 

PRINT 



VALUES 

VALLES 



4500 
5000 
5 5 00 



TIME DIFFERENCES 
CALCULATE SOURCE DEPTH, S-R RANGE 
PRINT DEPTH AND RANGE 
CALCULATE REFLECTION PATH DISTANC 



PS. 



PRINT DISTANCES 
YCU FINISHED? (Y OR N)— ” 



INPUT Z 4 $ 
PRINT "SAME 
INPUT Z$ 

IF Z $="N " 
Z2$="Y" 

GOTO 127 
REM 
REM 

IF Z $="N " 



6000 
"ARE 
Z$ 

»Y" THEN STOP 
"SAME CATA? (Y 



OR N) — 1 



PARAMETERS? 
GOTO 115 



(Y OR N) — " 



PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

PRINT 

INPUT 

PRINT 

PRINT 

PRINT 

INPUT 



*** INI TIALIZATION 
GOTO 525 



PROCEDURE* 5 ** 



"AUTOPEAK" 

"SAMPLING FREQUENCY MUST 
" GREATEST FREQUENCY OF 



BE FOUR TIMES 
INTEREST." 



THE" 



"NUMBER OF POINTS PER BLOCK — " 

N2 

"MINIMUM TIME DIFFERENCE BETWEEN DIRECT 
» AND FIRST REFLECTED PATH RECEPTION ( 
PRINT « SAMPLES) — " 

D8 



PATH" 

IN" 
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# -K- -Sf- -H- -if- 



535 PRINT "NUMBER CF PCINTS TO BE USED FROM EACH FL CCK" 
540 PRINT " (NUMBER CF POINTS PLUS DELAY MUST EE LESS" 
545 PRINT " THAN NUMBER CF POINTS IN BLOCK) — » 

55C INPUT N 1 

555 PRINT "MAXIMUM LAG — " 

560 INPUT LI 

562 PRINT "PRINT PEAK VALUES? (Y OR N)--" 

563 INPUT 23$ 

565 IF Z$="N" THEN PRINT "CONTINUING WITH C AL CU L AT I CN S" 

5 7 C IF Z $= "N" THEN RETURN 

575 PRINT "IS THIS THE 1ST BLOCK OF A MULTIPLE RUN? (Y/N) 
5 8 C INPUT Z2$ 

585 PRINT "SAMPLING FREQUENCY — » 

590 INPUT SI 

595 PRINT "SPEED OF SOUND (M/SEC) — " 

600 INPUT C 

605 PRINT "DEPTH OF WATER (M) — " 

61 C INPUT H 

615 PRINT "DEPTH OF RECEIVER t M ) — " 

620 INPUT R 
625 RETURN 
630 REM 

1000 REM ***ROUTINE TO READ DATA FROM TAPE*** 

1004 PRINT "SWITCH TO HIGH SPEED." 

1005 IF Z2$="Y" THEN INPUT ON (2) V$ 

1010 K 1 =0 

1015 FOR 1=1 TO (N2/8) 

1020 INPUT ON (2) V$ 

1025 FOR J=1 TO 48 STEP 6 
1030 Vl$=" " 

1032 V1S(1, 6)="" 

1035 FOR K= 0 TO 5 

1040 IF V$( J+K, J+K)=" " GOTO 1050 
1045 V1$=V1$+V$( J+K, J+K) 

105 C NEXT K 

1055 A ( K 1 ) = VA L ( V 1 $ ) 

1060 K 1 =K 1+ 1 
1065 NEXT J 
107 C NEXT I 

1030 PRINT "CONTINUE?" 

1035 INPUT Zl$ 

1090 PRINT "CONTINUING WITH CALCULATIONS." 

1095 RETURN 
1100 REM 

1300 REM ***8UILD OTHER ARRAY*** 

1305 FOR 1 = 0 TO Nl-1 
1310 B ( I ) =A ( I +D8 ) 

1315 NEXT I 
1320 RETURN 
1325 REM 

1500 REM ***CROSS-C GRRE LAT I ON ROUTINE*** 

1505 A 8 = 0 
1510 B 8 = 0 

1515 FOR 1=0 TO Nl-1 
1520 A8= A8+A ( I) 

1525 B 8= E8+B ( I) 

1530 NEXT I 
1535 AS=A8/N1 
1540 B8=B8/N1 
1545 FOR 1=0 TO LI 
1550 N9=N1-I 
1555 S8=C 

1560 FOR J = 0 TO N9-1 
1565 I 9 = J + I 

1570 S8=S8+( A ( J )-A8)*( BU9 )-B8) 

1575 NEXT J 
1580 Y ( I )=S8/N9 
1585 NEXT I 
1590 RETURN 
1 5 9 c REM 

2000 REM ***DET ERM I ME INTERMEDIATE PEAK VALUES*** 

2005 IF Y ( 1 ) > Y ( 0 ) THEN E=1 
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2010 
2015 
20 2 C 
2025 
2030 
2035 
204C 
2045 

205 C 
2055 

206 C 
2065 
2070 
2075 
208 C 
2085 
205 C 
2095 
2500 
2505 
25 1 C 
2515 
2520 
2525 
253C 
2535 
2540 
2545 
2550 
2555 
2560 
2565 
300 C 
3003 
3005 
3010 
3015 
3020 
3025 

3030 

3035 

3040 



3045 
305 C 
3C55 
3060 
3065 
3070 

3075 

3076 

3077 

3078 
3075 
3080 
3031 
30 6 2 
3085 
3050 
3055 
3100 
3105 
3110 
3115 
3120 
3125 
313C 
3135 
3137 
3140 
3145 



IF VI1XYI0I THEN E=-l 
P9=-l 

FOR 1=2 TO LI 
IF E=(-l ) GOTO 2 06 C 
IF Y { I ) > Y ( 1-1) GOTO 2085 
P5=P9+1 

P ( PS , 0 )= 1-1 +C8 
P <PS,1 ) = Y( I - 1 ) 

E = -E 

GOTO 2085 

IF YIIXY(I-l) GOTO 2085 
P 5=P9+ 1 

P <P9,0) = I-1 + D8 
P ( P 5 , 1 ) = Y ( 1-1) 

E=-E 
NEXT I 
RETURN 
REM 

REM ***PRINT INTERMEDIATE PEAK. VALUES*** 

PR INT 

PRINT "—PEAK VALUES — " 

PRINT 

FOR 1=0 TO PS 

IF 1/5= I NT( I /5) GOTO 2540 

PRINT P( I, 1) t 

GOTO 2550 

PRINT P ( 1,1) 

PRINT 
NEXT I 
PRINT 
RETURN 
REM 

REM ***DETERMINE TWO PEAK VALUES-** 

Z8 = G 
E = C 
K 1 =- 1 

FOR 1=0 TO PS-3 
IF E=1 GOTO 3040 

IF A6S(P ( I, 1 ) )>ABS< P( I+lt 1 ) ) OR ABS(P(I,1))> 

ABS ( P( 1+2,1 ) ) GOTO 3075 
IF ABS(P<I,1))> ABS(P( 1+3, 1) ) GOTO 3075 
E=1 

IF ABS(P(I,1 ))<ABS(P( 1 + 1,1) OR A BS t P ( I , 1 ) 1 < 

A BS ( P ( 1 + 2,1) ) GOTO 3075 
IF ABS(P( 1,1 ))<ABS(P( 1+3,1)) GOTO 2C75 
E=C 

K 1=K 1+ 1 

M ( K1 ,0 ) =P ( I ,C)/S1 
M(K1,1)=P(I,1) 

M ( K1 » 2 )= I 

NEXT K 

IF K1>=0 GOTO 3081 

PRINT "CHANGE LAG — MAXIMUM LAG= " 

INPUT LI 
Z 8= 1 
RETURN 

REM DETERMINE TWO HIGHEST PEAKS 
M8=M ( 0, 1 ) 

I 5 = C 

FOR 1=1, K1 

IF ABS ( M ( I , 1 ) )> AB S ( M8 ) THEN 19=1 
IF ABS(M(I,1))> A 3 S ( M 8 ) THEN M8=M(I,1) 

NEXT I 
R ( 0 , 0 ) = M 8 
R(C,1)=MC 19,0) 

R( C,2) = M( 19, 2) 

T g _ r 

IF 19=0 THEM 18=1 
M€ = M( I 8, 1) 

IF Kl= I GOTO 3165 
FOR 1=18+1 TO K1 
IF 1=19 GOTO 3160 
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31 5 C IF ABS ( M ( I , 1 ) )> ABS ( M8 ) THEN 18=1 
3155 IF A6S(M( I , 1) )> ABS ( M8 ) THEN M8=M(I,1) 

3 16 C NEXT I 

3165 R( 1,0)=M8 

31 7 C R(1,1)=M(I8,C) 

3175 R( 1,2)=M( 18, 2) 

3 1 8 C RETURN 
3185 REM 

4000 REM ***PRINT TIME DIFFERENCES*** 

4005 PRINT 
40 1 C PRINT 

4015 PRINT "TIME DIFFERENCE BETWEEN DIRECT AND SURFACE 
REFLECTED PATHS" 

402C PRINT " 

4025 PRINT TA8( 9 ) ;"PEAK VALUE" 

4027 A3 = R ( 0 , 1 )*10C0 

4030 PRINT USING » 33.53333 MSEC", A3 

4035 PRINT 
404C PRINT 

4045 PRINT "TIME DIFFERENCE BETWEEN DIRECT AND 6CTTCM 
REFLECTED PATHS" 

4050 PRINT " 



4055 PRINT TAB19) ;"PEAK VALUE" 

4057 A 3 = R ( 1 ,1 )*10C0 

406C PRINT USING » aa. 53333 MSEC", A3 

4065 PRINT 

40 7 C PRINT 

4075 RETURN 

4080 REM 

450C REM ***CALCULATE SOURCE DEPTH AND S-F RANGE*** 

4505 FOR 1=1 TO 1 

45 1 C G C = R ( 0 , I )* (C*R( 1 , I ) ) 2 

4515 G 1=C 2*R(0, I )*R( 1 , I)+4*H 2-4*H*R 

452C G2=R { 0, I )*G1 

45 25 G3=4*(R( 0, I ) * ( R-H ) -R*R (1,1) i 
4530 S ( I-l)=( G0-G2)/G3 

4535 C< I-l)=2*R*S(I-2 )/ CC*R (0, I ) )-C*R (C, I)/2 
454 C NEXT I 
4545 RETURN 
4550 REM 

500C REM. ***PPINT DEPTH AND DISTANCE VALUES*** 

5005 PRINT "RANGE OF SOURCE FROM RECEIVER AND SOURCE DEPTH 
IN METERS:" 

5010 PRINT " 



5015 PRINT TA B ( 9 ) ;"PEAK VALUES" 

5020 PRINT TAB( 5) ;"DEPTH"; TAB( 19) ; "RANGE" 

5030 PRINT USING " 33.53353 ",S(0),D(0) 

5035 PRINT 
504 C PRINT 
5045 PRINT 
505 C RETURN 
5055 REM 

5500 REM ***C ALCULAT E DISTANCES CF REFLECTED PATHS** 

5505 REM CALCULATE TRANSVERSE SOURCE-RECEIVER CISTANCE 
5510 T=SQR (D ( 0 ) 2-(R-S(C)) 2) 

5515 REM CALCULATE DISTANCE OF SURFACE REFLECTION PATH 
55 2 C X=S(0)*T/(R+S(0) ) 

5525 Y= F*T / ( R + S ( 0 ) ) 

5530 U=SQR( X 2+S(C) 2) 

5535 W= SQR ( Y 2+R 2) 

5540 DO=U+W 

5545 REM CALCULATE DISTANCE OF BOTTOM REFLECTION PATH 
5550 A=(H-S(0))*T/(2*H-R-S(C) ) 

5555 B=(H-R)*T/(2*H-R-S(0) ) 

5560 E = SQR ( A 2+(H-S(0)) 2) 

55o 5 F = SQR ( B 2+(H-R) 2) 

5570 D1=E+F 
5575 RETURN 
5580 REM 
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600 C 

6005 

6010 

6015 

6020 

6025 

603 C 

6035 

6040 

6045 

605 C 

6055 

606 C 

6065 



REM ***PRINT PATh DISTANCES*** 

PRINT "SURFACE REFLECTION PATH DISTANCE IN METERS:" 

PRINT " " 

PRINT USING " a/Sa. 235)0)2" , CO 

PRINT 

PRINT 

PRINT "BOTTOM REFLECTION PATH DISTANCE IN METERS:" 
PRINT « " 

print using « aaa.aaaaa ", di 

PRINT 

PRINT 

RETURN 

REM 

END 
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COMPUTER PROGRAM DEVERB 



* DEVERB 10/21/77 * 

* JEANIE SAVAGE, PROGRAMMER * 

* SPECIFICATIONS BY RICK BOSTIAN * 

* IN THIS PROGRAM DEVER3E RAT I ON * 

* IS PERFORMED IN THE FRE- * 

* QUENCY DOMAIN. * 

* LAST CORRECTION: 12/06/77 * 

£ # # * $ * £ =|£ * if ^ * * if. ^ * 5k aj: # jj; ^ £ £ £ £ £ £ $ ^ ^ jjt 2js * 

INTEGER XGRID 

I NT EGER *2 IZ2,IZ3,IZ4,IZ5,IZ6,YES,IZ7 
DATA YES/ 1 Y • / 

DIMENSION I STACK (20) , A ( 10 00 ) , 3 ( 1 OCO ) , I ARY ( 1C00) 
DIMENSION X(IOOO) , Y( 1000) 

DIMENSION FINAL ( 1 COO, 2) , I BEG( 20 ) , I END ( 2 0 ) , I CASE ( 20 ) 
COMMON SF,ISF, IBM, THETA, SIGN A, GAMMA, C,D,RC0EFF,N2, 
SISIZE,DB,DS,NBLK, IZ3, I Z 4 

T *V "T* T T* 'T* 'I' T *T* 'l* “1' 'T' T “T* V V ■'f* I 'T "l* 

INITIALIZATION ROUTINE 

I READ=0 
PRINT 500 

500 FORMAT! 'O' , 'DEVERB' ) 

PRINT 510 

510 FORMAT (• O' , 'NUMBER OF POINTS PER SIGNAL (PCWER OF 2)'. 
S' (15) — •) ' 

READ 520, N2 

520 FORMAT (15) 

PRINT 530 

530 FORMAT ( ' ','IS THIS THE FIRST SIGNAL OF A MULTIPLE', 

S' RUN? (Y/N)-' ) 

READ 540, I Z 2 

540 FORMAT ( A 1 ) 

PRINT 550 

550 FORMAT! ' ', 'SAMPLING FREQUENCY (F9.3) — •) 

READ 560, SF 

560 FORM AT ( F9 . 3 ) 

PRINT 570 

570 FORMAT ( ' ', 'DIRECT PATH DISTANCE (F9.5) — ') 

READ 580, D 

580 FORM AT ( F9 . 5 ) 

PRINT 590 

590 FORMAT ( ' ', 'SURFACE PATH DISTANCE IN METERS ( F9 . 5 ) — ') 
READ 580, DS 
PRINT 600 

600 FORMAT ( * ', 'BOTTOM PATH DISTANCE IN METERS (F9.5) — ') 
READ 580, DB 
PRINT 610 

610 FORMAT ( ' ', 'SURFACE REFLECTION TIME IN MSEC (F9.5 ) — •) 

READ 580, TS 
TS=TS/1 000 
PRINT 620 

620 FORMAT ( ' ', 'BOTTOM REFLECTION TIME IN MSEC (F9.5)— ') 
READ 580, TB 
TB=TB/1 000 
PRINT 630 

630 FORMAT ( ' ', 'BOTTOM REFLECTION COEFFICIENT (FS.5) — ') 

READ 58 C , RCOEFF 
WR I TE ( 6 , 71 1 ) 

711 FORMAT ( ' ','RMS WAVE HEIGHT (F9.5) — •) 

READ 580, SIGMA 
PRINT 720 

720 FORMATt' ', 'SURFACE ANGLE IF INCICENCE (IN RADIANS)', 
S' (F9.5) — ' ) 

READ 580, THETA 
PRINT 730 

730 FORMAT ( ' ', 'BOTTOM PHASE SHIFT (F9.5) — •) 

READ 580, GAMMA 
PRINT 640 
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640 FORMAT ( ' ', 'SPEED OF SOUND (F9.3) — ') 

READ 56 Ot C 
PRINT 650 

650 FORM AT ( * ', 'FREQUENCY PLOT? (Y/N)--') 

READ 540, I Z 3 
PRINT 655 

655 FORMAT ( ' ','TIME PLOT BEFORE FFT? (Y/N) — •) 

READ 540, I Z 7 
I F ( IZ7.NE .YES) GOTO 659 
PRINT 657 

657 FORMAT ( * * ,' NUMBER OF POINTS TC BE PLOTTED (15) — •) 
READ 520, INUM1 

659 PRINT 660 

660 FORMAT ( ' ','TIME PLOT AFTER FFT? (Y/N) — •) 

READ 540, I Z 4 

I F ( IZ4.NE.YES) GOTO 672 
PRINT 657 
READ 520, INUM2 
672 PRINT 675 

675 FORMAT ( ' ' , * ALL BLOCKS? (Y/N) — •) 

READ 540, I Z 5 

IF( IZ5. EQ.YES) GOTO 1000 

KPTR=0 

PRINT 680 

680 FORMAT ( ' ', 'SPECIFIC BLOCKS (12) (INPUT 9<= WHEN', 

S' FINISHED)—' ) 

690 READ 7 00 , I TE MP 
700 FORMAT (12) 

IFUTEMP.EQ.99) GOTO 1000 
KPTR=KPTR+1 
I STACK( KPTR)=ITEMP 
GOTO 69 C 

READ DATA TAPE 

4 * 

1000 I F ( IREAD.EQ.l ) GOTO 2000 
PRINT 1005 

1005 FORMAT ( * ', 'READY TO READ DATA TAPE') 

I F ( IZ2.NE.YES) GOTO 1020 
READ (5,1010 KEND 
101C FORMAT (816) 

1020 DO 1030 1=1, N2, 8 
I TEMP= 1+7 

READ( 5, 1010 ) ( IARY(J) ,J = I ,ITEMP) 

1030 CONTINUE 

* * * * £ $ # * * * * * * 4= -if * * * * * $ * * * * * * * £ 

DETERMINE DIVISION BOUNDARIES 

iytf. If. if. -3y% 

PRINT 1501 

1501 FORMAT ( ' ' , « CONT INU I NG WITH CALCULATIONS') 

15 OC DO 1510 1 = 1, N2 
T I ME= ( I -1 >/SF 
IF(TIME.LT.TS) GOTO 1510 
I SF= I 
GOTO 1520 
1510 CONTINUE 
I SF=N2+ 1 
GOTO 1540 

1520 DO 1530 1=1, N2 
T I ME= ( 1-1 ) / SF 
IF(TIME.LT.TB) GOTO 1530 
IBM=I 

GOTO 1550 
1530 CONTINUE 
1540 IBM=N2+1 
C 

C ***DETERMINE NUMBER OF POINTS IN EACH SECTOR*** 

1550 NPT1= I SF-1 
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NPT2 = IBM- ISF 

I F ( IBM. IT. I SF) NPT2= I SF- I BM 
NPT3=N2-I BM + 1 

I F( IBM. LT . I SF) NPT3=N2-ISF+1 

Ota?*. £ ***:£■* 

DETERMINE BLOCK SIZE 

JU U# Uf U# U# 4# 4o 4L» <JU 4# 4# 4# X X X X X 

-*r ^ *4* ^ T 1 ■'T- -r* ^ *r t t 

***FIND SMALLEST AND MIDDLE NUMBER OF POINTS IN SECTOR 
1700 ISMALL=MINC(NPT1 ,NPT2) 

ISMALL=MIN0(NPT3, I SMALL) 

I F { ISMALL.EC.NPT1 ) MI DDLE=M IN0< NPT 2, NPT 3) 

IF( ISMALL.EC.NPT2 ) M I DDL E=M INO ( N PT 1 , NPT 3 ) 

I F ( I SMALL .EC.NPT3 ) M I CDLE=M INO ( N PT 1 , NPT 2 ) 

***DETERMINE NUMBER OF POINTS PER BLOCK*** 

2000 I F ( IABS (MIDDLE/2-ISMALL) .GT . { M I DOLE- 1 SMALL ) ) 

£1 SIZE=ISMALL 

IF( ISIZE.EQ.ISMALL) GOTO 2200 
DO 2010 K=l,30 

I F ( ISMALL.EC.MIDDLE/K ) GOTO 2020 

I F( ( I SMALL. EQ.MI DDLE/K) . LT. ( I SMALL-MI DDL E/ < K+l) ) ) 

SGCTO 2020 
2010 CONTINUE 
2020 I S I Z E= M I DDLE/K 

***DETERMINE NUMBER OF BLOCKS PER SECTOR*** 

2200 NBLK1=NPT1/ISI ZE 
NBLK2=NPT2/ISI ZE 
N8LK3=NPT3/ISI ZE 

***DETERMINE NUMBER OF POINTS SKIPPED EACH SECTOR*** 
2300 NSKI P1=NPTI-NBLK1*ISIZE 
NSKIP2=NPT2-NBLK2*I SI ZE 
NSKIP3=NPT3-MBLK3*ISI ZE 

***PR I N T SECTOR INFORMATION*** 

WRITE(6 ,2405) ISIZE 
2405 FORMAT ( *0' ,' ISI ZE=' , 13) 

I T 1 = 1 
I T2 = 2 
I T 3 = 3 

WR ITE ( 6 , 241 C ) I Tl , NBLK1 ,NSK IP1 
2410 FORMAT { '0' , 'SECTOR ',11,' CONTAINS ',12,' BLOCKS,', 

$ 2X , 13 , ' POINTS SKIPPED') 

WRITE(6,241C) I T2 , N0 LK2 , NSK I P 2 
ViRITE(6,241C) I T3 , NBL K3 , NSK I P3 

=* $:{: ^ * £ # 

BUILD STACK IF PROCESSING ALL BLOCKS 

*•!{;* £:$: ^ ^ ## £ * * £ * $ * 5 ?-- # 3 : #* # 

500 IF( IZ5.NE.YES) GOTO 2700 
NUMBLK=N2/I SIZE 
DO 2510 1=1 ,NUMBLK 
ISTACKI I ) = I 
510 CONTINUE 

KPTR=NUMBLK 

*** ** £ * £ * £ 

DETERMINE BEG £ END POINTS £ CASE 

700 16=0 

***BLOCKS IN 1ST SECTOR*** 

ITEMP=NSKIP1+1 
DO 2710 1=1 ,NBLK1 
16=16+1 

I BEG( 16 )= I TEMP 
I END( I6)=ITEMP+ISIZE-1 
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I CASE ! I 6 ) =1 
ITEMP= I END( I 6 ) + 1 
2710 CONTINUE 
C 

C ***BLOCKS IN 2ND SECTOR*** 

I TEMP= I SF 
I TEMP1 = NS KI P2 
DO 2720 1=1 i NB LK2 
I 6=16+1 

I B EG ( 16 )= ITEMP 
I END( 16 )= ITEMP+I S I ZE-1 
I T EMP=I END( 16 ) + l 
I CASE ( I 6 ) = 2 

I F ( ISF.GT.IBM) ICASE!I6)=3 
I F ( ITEM P2 . EQ.O ) GOTO 2720 
I TEMP=I TEMP+1 
ITEMP1=ITEMP1-1 
272C CONTINUE 
C 

C ***BLOCKS IN 3RD SECTOR*** 

DO 2730 I=1,NBLK3 
16=16+1 

I B EG ( 16 )= IBK + ! 1-1 )*ISIZE 
IF{ ISF.GT .IBM) IBEG( I 6 ) = I SF + ( I-1)*ISIZE 
I END( 16 ) = I BEG( I6J+ISIZE-1 
ICASE ( I 6 ) =4 
2730 CONTINUE 
C 

C V< Vr *4 "A* U* U/ 4* *V 4( O, OL> JU 4/ 4/ 4# 4> 

'r'rT'r^'T'P'r^T -r* 4 s -V 1 -¥■ -r- -r* 

C PRINT TIME PLOT BEFORE FFT 

C 

2900 I F ( IZ7.NE.YES) GOTO 3000 

IF! INUM1.GT .1000) GOTO 2910 
NUM=I NU Ml 
K P LOT S= 1 
GOTO 2920 

2910 KPLOTS= INUM1/100Q+1 
NUM=1000 

2920 DC 2940 I=1,KPL0TS 

I F ( ( I .EQ.K PLOTS ) .AND. (I .NE. 1) ) NUM=INUM1- (K PLOTS- 1 ) 
6*1000 

DC 2930 J=1,NUM 

Y ! J ) =FL 0 AT! I ARY { ( I - 1 ) *1 000+ J ) ) 

X<J)=<( 1-1 )*1000+J-1 ) /SF 



2930 


CONTINUE 
XGRID=NUM/1 C 
I F( MOD! NUM, 10) .ME .0) 


XGRI D=XGR I C+ 1 




CALL PL OTDV ( X, Y , NUM, 


X GRID, 3 ,NBLK ) 


2940 

c 


CONTINUE 




c 


U# JLfJU JU JU V** **^ 

^ ^ oA rp 




C 


PROCESS BLOCKS 




C 


yp op T o' •nr ^ ^ or o' 





C 

300 C DO 3110 1=1 ,KPTP 
NBLK= 1ST ACK ( I ) 



C 

C ***DETERM INE BOUNDARIES AND CASE CF BLOCK*** 

I 8EG0=I BEG ( NBL K ) 

I END0=I END ( NBLK ) 

I C ASEO= ICASE(NBLK) 

C 

C ***PROCESS BLOCK*** 

PRINT 3015, IBEGO, IENDO, ICASEO 
3015 FORMAT!///* * , « I B EG= * , 14, 4X , * I END= * , 14, 4X , • I C ASE= • , 1 1 ) 
K = 1 

DO 3020 J= I BEGO , I ENDO 
A(K)=FLOAT( IARY(J) ) 

B ( K )= 0. 0 
K = K+1 

3020 CONTINUE 
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CALL SIGNAL ( A , B , F INAL , I BE GO , I C A S EC , I S I Z E J 
IF(< I .NE.l) .OR. (NSKIP1.EQ.O) ) GOTC 3040 
OC 3030 J=1,NSKIP1 

$LOPE=F INAL ( I BE G ( I J+l ,1 ) -FINAL! I BEG (I ) , 1 ) 
FINAL(J,1)=-SL0PE+FINAL(I8EG( I) ,1) 

SLOPE=F INAL <IBEG(I)+1,2)-FINAL(IBEG(I),2) 

F INAL( J ,2 )=-SLOPE+FIN AL ( IBEG( I ) , 2 ) 

3030 CONTINUE 

3040 I F ( (I .GT.NBLK1+NBLK2) .OR. (I .LE.NBLK1+1) ) GOTO 3110 
IF(NSKI P2.EC.0) GOTO 3110 
DO 3050 J=1,NSKIP2 
I TEMP = I END ( NBLK1 + I ) + l 

TEMP=F I NAL ( ITEMP+1,1 J-FINAL (ITEM P-1,1) 

FINAL (I TEMP ,1)=FINAL( ITEMP- 1 , 1 ) +T EMP 
TEMP=F I NAL ( ITEMP + 1 . 2 ) -F INAL < IT EMP-1 , 2 ) 
FINAL(ITEMP,2)=FINAL( ITEMP- 1 , 2 ) +T E MP 
3050 CONTINUE 
3110 CONTINUE 
C 

C -JU >u XU XU XU XU «JU XU XU XU XU XU J* -at Xu XU «l* J. 

C PRINT TIME PLOT AFTER FFT 

C «V XU XU "Jr XU X X -»V» XU XU Uf XU X X X X X "Jr X 

t* ^ -v* **■> or -r- or -r -t* -t 

C 

3300 I F ( IZ4.NE .YES) GOTO 3500 
DO 3308 1=1 ,KPTR 
NBLK= I STACK ( I ) 

N= I S I ZE 

IBEGO=IBEG(NBLK J 
I END0=I END ( N8LK ) 

K- 1 

DC 3304 J=I BEGO , I ENDO 
A ( K ) =F I NAL ( J, 1 ) 

8<K)=FI NAL( J,2) 

K =K + 1 

3304 CONTINUE 

CALL CFFT2< A,B,N ,N,N, 1 ) 

DO 3306 J-l , IS I ZE 
F IN AL ( I BEG0 + J-1 , 1 ) = A( J) 

3306 CONTINUE 
3308 CONTINUE 

IF(NSKIPl.EC.O) GOTO 3313 
DC 3312 1 = 1 ,NSK I PI 

FINAL (I ,1)=(FINAL(IBEG(1),1 >/(NSKIPl+l) )*I 
3312 CONTINUE 

3313 I F { NSKI P2 .EG.O ) GOTO 3315 
DO 3314 1=1, NSK I P2 
I TEMP=( NBLK1+I ) +1 

T EMP=F I NAL ( IT EMP + 1 ,1 ) -FI NAL (I T EM P-1, 1 ) 

FINAL! I TEMP, 1 )=FINAL< I TEMP-1 , 1 ) +T EMP 
3314 CONTINUE 

3315 I F ( INUM2 . GT . I END ( KPTR) J I NUM 2= I END (KPTR) 
IFUNUM2.GT.1000) GOTO 3310 
NUM=I NUM2 
K PLOTS= 1 
GOTO 3320 

3310 KPLOTS= IMUM2/1 000+1 
NUM=1 00 0 

3320 DO 3340 I=1,KPL0TS 

I F ( { I . EQ . KPLOT S ) .AND. (I. NE.l)) NUM=I NUM 2- ( K PLOTS- 1 ) 
8*1000 

DC 3330 J =1 , NUM 
Y ( J ) =F I NAL ( ( I -1 ) * 1 000+J , 1 ) 

Y ( J ) = Y( J ) /FLOAT ( I SI ZE) 

X ( J)=( ( 1-1 ) *1000+J-1 ) /SF 
3330 CONTINUE 

XGRID=NUM/10 

I F ( MOD ( NUM,10).NE.0) XGRI D=XGRI 0 + 1 
CALL PLOTDV(X, Y,NUM,XGRID ,2 ,NBLK ) 

3340 CONTINUE 
C 

C * x. xux — u xu „u 

C CONCLUSION 
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c 

c 

3500 

3510 



C 

C 

C 

C 



C 

C 

c 

c 

c 

100 

c 

c 

c 

c 

c 

500 



c 

c 



c 

c 

510 



520 

C 

C 

C 

C 

C 

1000 

C 

C 



-x, a^> -jUvi. j. -ju ju 

t* *r -y y *v ~r- 

PRINT 3510 

FORMAT ( ' 1 , 'ARE YOU FINISHED? ( Y/N ) — * ) 

READ 540, I Z 6 

I F ( IZ6. EQ .YES) STOP 

I READ=1 

GOTO 672 

DE3UG SUBCHK 

END 



****SIGNAL SUBROUTINE**** 

SLBROUT INE SIGN AL < A , B , F I NAL , I B E GO , IC AS E 0 , N ) 
DIMENSION F INAL ( 1 000 * 2 ) 

DIMENSION X( 1000) t Y( 1000) ,A( 1000) ,B( 1000) 

INTEGER*2 YES,IZ3,IZ4 

COMMON SF , ISF, I BM,TH-TA, SIGMA,GAMKA,C, D,RC0EP C ,N2, 
SI SI ZE » DBt OS ,N8LK, IZ3,IZ4 
INTEGER XGRID 
REAL KO 

DATA YES/’Y ’/ 



PERFORM FFT 



CALL CFFT2(A,B,N,N,N,-1) 

-V **■» ^ 4 » «->- J- -JU J» J/ JU «JL. 4 > 
TTT'ri*T'i''T''(‘Ti''rT / i' -v -t* *0 

PERFORM CORRECTIONS 



DO 520 1=1, N 

FINAL! I 8EG0+I-1 , 1 ) =A ( I ) 

FINAL! I 3EG0 + I-1 ,2 )=B( I) 

FREQ= ( I -1 ) * SF/N 

PI=3. 14159 

K0=2.*P I*FREQ/C 

IF! ICASEO.EQ.l ) GOTO 520 

IF! ICASEO.EC-3) GOTO 510 

***CORRECT I ON FOR SURFACE REFLECTION*** 

G=( !4.*PI*SIGMA*FPEQ/C)*C0S(THETA ) )**2 
I TEMP= ( IBEG0+I-1)-(ISF-1) 

TMAG=SQRT (F INAL ( I TEMP 1 1 ) **2 + FI NAL ( ITEMP , 2 )**2 ) 
TPHAS6=ATAN2(FINAL( ITEMP, 2) , F IN AL ( ITEMP , 1 ) ) 
S=D/DS*EXP( -G/2. ) *TMAG 
SPHASE=TPHASE-! DS-D)*KO-PI 

FINAL! I BEGO + I-1 , 1 )=FINAL( I BEGO+ I - 1 , 1 ) -S *COS (SPHASE) 
FINAL! I BEGO + I-1 , 2 ) = FINAL( 1 6 EGO + 1 - 1 , 2 ) -S*S I N ( SPH AS E ) 
IF! ICAS EO . EC .2 ) GOTO 520 

***CORR ECT I ON FOR BOTTOM REFLECTION*** 

I TEMP= ( IBEGC+I-1 ) — ( I 8 M — 1 } 

TMAG=SQRT INAL (ITEMP, 1 ) **2 + F I N AL ( I T EMp , 2 )**2 ) 
TPHASE= AT AN2! FINAL! ITEMP, 2) , F I N A L ( IT EMP , 1 ) ) 
S=RCOEFF*D/DB*TMAG 
SPHASE=TPHASE-(DB-D)*KO+GAMMA 

FINAL! I BEGO+I-1, 1 )=T INAL ( 13 EGO+ I- 1 , 1 ) -S*C OS ( SPHAS E ) 
F INAL! I BEGO + I-1, 2 )=F INAL! IBEGO+ I - 1,2 )-S*S IN (SPHASE) 
CONTINUE 

* * *** *** * * * **** ***'**:$: * *** * 

PRINT FREQUENCY AND TIME PLOTS 

IF! IZ3.NE .YES) GOTO 1020 

***FREQUENCY PLOT*** 

K = 0 
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I TEMP=N /2 
DO 1010 I=1,ITEMP 
K=K + 1 

Y(K)=FINAL ( IBEGO+I-1 , 1) **2+ FINAL ( I8EG0+ 1-1 , 2) **2 
Y(K)=10*ALOG10( Y( K)*D) 

X ( K)=( K-l ) *SF/N 
1010 CONTINUE 
NUM=K 

XGRID=NUM/1 C 

I F ( MOD ( NUM, 10) .NE.O) XGRI D=XGR I C+ 1 
CALL PL0TDV(X,Y,NUM,XGRID,1 ,NBLK) 

1020 RETURN 

CEBUG SU6CHK 
END 
C 
C 
C 

C ****PI_OT SUBROUTINE**** 

C 

SUBROUTINE FLOTOV ( X, Y, N, XGR I D, M t Ne ) 

INTEGER D,XGRID, YGRID,4XIS 

DIMENSION Y (1000) ,C7(101),0(6),X(1000),KAXIS(51> 
DATA IDASH/1H-/, ISTAR/1H*/, ID0T/1H./ 

DATA IBAR/1HI/ , IPLUS/1H+/, I BLANK/lh /,IX/1HX/ 

AXI S=51 
YGRID=6 
XGR I D=X GR I D+l 
212G N1=N-1 

Y 6 = Y ( 1 ) 

Y 1=Y ( 1 ) 

DO 2200 1=1 ,N 

I F(Y6-Y( I ) .GE.0.0 ) GOTO 2180 

Y 6=Y ( I ) 

2180 I F ( Yl-Y( I ).IE.0.0) GOTO 2200 
Y 1=Y ( I ) ' 

2200 CONTINUE 

S=Y1* (AXIS-1)/ (Y6-Y1 ) 

X1=X( 1) 

XI 0=X( N ) 

0( 1)=Y1 
0 ( 6 ) = Y6 
I I X=XGR ID-1 
DO 2410 1 = 1, 1 1 X 
C7( I)=X( ( I-l)*10+l) 

2410 CONTINUE 

C7(XGRID)=C7(XGRID-1 ) -»1 0* ( X ( 2 ) - X ( 1 ) ) 

I FCN.EQ. (XGR ID- 1 )*10) C7(XGRID)=X(N) 

I IY=YGR ID-1 
DO 2440 1=2, IIY 

0 ( I ) = ( F LO AT ( 1-1 ) * ( Y6 — Y1 ) / FL 0 AT ( YGR I D- 1 ) )+Yl 
2440 CONTINUE 

WRITE( 6,2460) 

2460 FORMAT ( ///, • • ) 

IF(M.NE.l) GOTO 2485 
PRINT 2470, NB 

2470 FORMAT ( *0 , ,32X, , BL0CK',1X,I2) 

2485 IF(M.EQ.l) PRINT 2486 

2486 F CRM AT ( * • t 27X,*DB ,, S VS. FREQUENCY') 

IF(M.EQ.2) PRINT 2488 

2488 FORMAT ( 1 ', 23X, ' VOLTAGE VS. TIME (AFTER FFT ) ' ) 

IF(M.EQ.B) PRINT 2487 

2487 FORMAT ( * ' , 23X , • V OLT A GE VS. TIME (BEFORE FFT ) ’ ) 

WRITE (6, 2500) ( 0 ( I ) , I = 1 , YGR I D ) 

2500 FQRMAT(9X»11(1PE10.2) ) 

S1=(X10-X1) /10.0*(XGRID-1) 

D= 1 
L 1 = 1 
L = 1 

1 2= I F I X ( -S+ 1 . 5 ) 

ITEMP=( XGRID-1 )*10+1 
DO 2900 I1=1,ITEMP 
IFCN.LT.il) GOTO 2510 
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YTEMP=( Y( ID *FLOAT<YGRID- 1 1 *10. 0/(Y6-Yl))-S 
I Y=I F I X I YTFMP+ 1 .5) 

2510 IFILl.GT.L) GOTO 2760 
DO 2650 IP=1,AXIS 
KAXISI IP) = IDASH 
2650 CGNTINUE 

DO 2680 1 = 1 , AXIS, 10 
KAXISI I )= IP LUS 
268 C CONTINUE 

IF(N.LT.Il) GOTO 2720 

I FI (Yl.LE .0.0) .AND.( 0.0.LE.Y6) ) KAXISI IZ ) = I COT 
KAXISI I Y ) = I STAR 

2720 WRITEI6 ,2725) C 7 ( D), ( KAXISI J) ,J = 1 ,AXIS) 

2725 FORMAT! 1PE13- 2 ,2X,115A1) 

L1=L1+1 0 
D = D + 1 
GOTO 2870 

2760 DO 2780 IP=1,AXIS 
KAXISI IP )=I BLANK 
2780 CONTINUE 

DC 2810 1=1, AXIS, 10 
KAXISI I )= 1 8 AR 
2810 CONTINUE 

IFIN.LT.il) GOTO 2860 

IF! (Yl.LE. 0.0) .AND. (0.0.LE.Y6) ) KAXISI I Z )=ICOT 
KAXISI IY)=I STAR 

2860 WRITE! 6 ,2 865 ) I KAXIS I J ) , J = 1 , AX IS ) 

2865 FORMAT! 15X, 115A1 ) 

287C L=L+1 

IFIL.GT .( XGRID-1 ) *10 + 2) GOTO 2910 
2900 CONTINUE 
2910 RETURN 

CE8UG SU3CHK 
END 
C 
C 

C ***FCURIER TRANSFORM SUBROUTINE*-* 

C 

SUBROUTINE CFFT 2 I A, B , NTOT , N , NSP AN, ISN ) 
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COMPUTER PROGRAM TDEVERB 



* $ ^ ifr if. * ^ ^ 

* TDEVERB 12/07/77 * 

* JEANIE SAVAGE, PROGRAMMER * 

* SPECIFICATIONS BY RICK BOSTIAN * 

* IN THIS PROGRAM D EVERBE RAT I ON * 

* IS PERFORMED IN THE TIME * 

* DOMAIN. * 

* LAST CORRECTION: 12/08/77 * 

£ ^ ****3r# =r ^ ^ ^ 5* !js 

INTEGER XGRID 

INTEGERS IZ2,IZ3,IZ4,IZ5,IZ6,YES,IZ7 
DATA YES/'Y •/ 

DIMENSION I STACK (20) ,A(1000),S< 10 00 ) , I ARY ( 1 COO ) 
DIMENSION X(IOOO) ,Y( 1000) 

DIMENSION FINAL ( 1000) ,IBEG< 20 ) , I END ( 20 ) , I C A SE ( 20 ) 

■JU Vf <JU d* *JU -JU -JU <v 'A* >Y "4. YY *L J. Vf 4. 4, 4. W 4/ 4, -JU» 

-i» ■nr T^^n'vi‘’ , T' , r*4r ;: *'r 7 r *r t 'r *r n* *r‘ m 

INITIALIZATION ROUTINE 

I READ=0 
PRINT 500 

500 FORMAT ( 'O' , 'DEVERB' ) 

PRINT 510 

510 FORMAT( 'O' , 'NUMBER OF POINTS PER SIGNAL (PCWER 0<= 2)', 
S' (15)—') ' 

READ 520, N 2 

520 FORMAT (15) 

540 FORMAT (Al) 

PRINT 550 

550 FORMAT (' 'SAMPLING FREQUENCY (F9.3) — •) 

READ 560, SF 

560 FORMAT ( F9 . 3 ) 

PRINT 570 

570 FORMAT (• ', 'DIRECT PATH DISTANCE (F9.5) — •) 

READ 580, D 

580 FORMAT ( F9. 5 ) 

PRINT 590 . 

590 FORMAT (• ', 'SURFACE PATH DISTANCE IN METERS <F9.5)--'j 
READ 580, DS 
PRINT 600 

600 FORMAT (• ', 'BOTTOM PATH DISTANCE IN METERS (F9.5) — •) 
READ 580, DB 
PRINT 610 

610 FORMAT (' ', 'SURFACE REFLECTION TIME IN MSEC (F9.5)--') 
READ 580, TS 
T S=TS/1 000 
PRINT 620 

620 FORMAT (• ', 'BOTTOM REFLECTION TIME IN MSEC (F9.5) — ') 
READ 580, TB 
TB=TB/1 000 
PRINT 630 

630 FORMAT (' ', 'BOTTOM REFLECTION COEFFICIENT (F9.5) — •) 

READ 580, RCOEFF 
WRITE(6 ,711 ) 

711 FORMAT ( ' ','RMS VnAVE HEIGHT (F9.5) — ') 

READ 580, SIGMA 
PRINT 720 

720 FORMAT ( ' ', 'SURFACE ANGLE OF INCIDENCE (IN RADIANS)', 
S' ( F9 . 5 ) — ' ) 

READ 580, THETA 
PRINT 730 

730 FORMAT ( ' ’, 'BOTTOM PHASE SHIFT (F9.5) — ') 

READ 530, GAMMA 
PRINT 640 

640 FORMAT (' ','SPEEC OF SOUND (F9.3) — ') 

READ 560, C 
PRINT 645 

645 F ORMAT ( ' ', 'BLOCK SIZE (15) — ') 

READ 520, ISIZE 
PRINT 650 
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650 FORMATl' ' , 'FREQUENCY PLOT BEFORE CORRECTIONS?', 

a* (y/n ) — ' ) 

RE AO 540, I Z 3 
PRINT 652 

652 FORMAT { * ’ , ' FREQU ENCY PLOT AFTER CORRECTIONS? (Y/N) — •) 
READ 540, I Z 2 
PRINT 655 

655 FORMAT ( ' ','TIME PLOT BEFORE CCRRECTIONS? (Y/N) — ’) 
READ 540, 111 
I F ( IZ7.NE.YES) GOTO 659 
PRINT 657 

657 FORMAT ( * », 'NUMBER OF POINTS TO BE PLOTTED (15) — ») 
READ 520, INUM1 
655 PRINT 660 

660 FORMAT ( ' ','TIME PLOT AFTER CORRECTIONS? (Y/N) — *) 
READ 540, I Z 4 
I F ( IZ4.NE.YES) GOTO 672 
PRINT 657 
READ 520, INUM2 
672 FRINT 675 

675 FORMAT (' ','ALL BLOCKS? (Y/N) — *) 

READ 540, I Z 5 

IF( IZ5.EQ.YES) GOTO 1000 

KPTR=0 

PRINT 680 

680 FORMAT ( ' ', 'SPECIFIC BLOCKS (12) (INPUT 99 WEEN FINIS 

690 READ 7 00 , I TEMP 
700 FORMAT (12) 

IF( ITEMP.EQ.99 ) GOTO 1000 
KPT R=KP TR+1 
I STAC K ( KPTR ) = ITEMP 
GOTO 690 

************** 

READ DATA TAPE 

*p 4*P -p -p p *p «Tp *p p p 

1000 IF ( IREAD.EG.l ) GOTO 2500 
PRINT 1005 

1005 FORMAT ( ' ', 'READY TO READ DATA TAPE') 

1010 FORMAT (816) 

1020 DO 1030 1=1, N2, 8 
I TEMP=I +7 

READ(5, 1010) ( IARY( J) ,J = I ,ITEMP) 

1030 CONTINUE 

DETERMINE DIVISION BOUNDARIES 
PRINT 1501 

1501 FORMAT ( ' • , • CONT I NU I NG WITH CALCULATIONS') 

1500 DO 1510 1=1, N2 
T I ME= ( I -1 ) / SF 
IF(TIME.LT.TS) GOTO 1510 
I S F = I 
GOTO 1520 
1510 CONTINUE 
I SF=N2+ 1 
GOTO 1540 

1520 CO 1530 1=1, N2 
, TIME=( 1-1 )/SF 

IF(TIME.LT.TB) GOTO 1530 
I BM= I 
GOTO 2500 
1530 CONTINUE 
1540 IBM=N2+1 



£ SjS >>5 # * £ £ * t- * $ $ £ * * * 5* # # * # * $ * # $ 3 * ^ * * £ 5ft * 

BUILD STACK IF PROCESSING ALL BLOCKS 
*****%*************** X************** 
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c 

250C 



25 1 C 

C 

C 

C 

C 

c 

2900 

C 

C 



I F ( IZ5.NE.YES) GOTO 2900 

NUKBLK=N2/I SIZE 

DO 2510 1=1 tNUMBLK 

I ST ACK ( I ) = I 

CONTINUE 

KPTR=NUMBLK 

■O* ^ ^ «JU%U O* «JU ^ J# «i 

T ‘T O' ^ ^ ^ 4 



- — — ' ' — — TT ’JU U# *J/ 

h nr -t> -nr -r -r 'r^r -tttt ^ 'i* *i+ 



PRINT TIME F-L0T BEFORE CORRECTIONS 






> «L U. U< V 

■ -nr " 



■ «•»* nr -nr* nr 






IF( IZ7.NE.YES) GOTO 3000 
I F ( IZ5.NE.YES) GOTO 2950 



***TIME PLOT FOP ENTIRE SIGNAL*** 

IF ( INUM 1.GT.1000) GOTO 2910 
NUM=I NU Ml 
KPL0T S= 1 
GOTO 2920 

2910 K PLOT S= I MUM 1/10 00 + 1 
NUM=1000 

2920 DC 2940 I=1,KPL0TS 

IF! ( I .EQ.KPL0TS ) . AND. (I .NE. 1) ) NUM-INUM1-! KFLCTS-1 ) * 
& 1 COO 

DC 2930 J=1,NUM 

Y(J)=FL0AT( I ARY ( ( 1-1) *1000+J) ) 

X(J)=( ( I-1)*1000+J-1)/SF 
2930 CONTINUE 

XGRI D=NUM/1 0 

I F ( M0D( NUMt 10J.NE.0) XGR I D=XGR I C+ 1 
N8LK=- I 

CALL PL0TDV!X,Y,NUM,XGRID,3,N3LK ) 

2940 CONTINUE 
GCTO 3000 



C 

c 



2950 



296C 



2970 

C 

C 

C 

C 

C 

C 

C 

3000 



* **T I ME PLOT FOR INDIVIDUAL BLOCKS**- 
DC 2970 I=1,KP*R 
NSLK= I STACK ( I ) 

DO 2960 J=1,ISIZE 

Y ( J ) = FL0AT ( I ARY ( N BLK- 1 ) *ISI ZE + J ) 

X( J)=( ( NBLK-1 )*ISIZE+J-1 )/SF 

CONTINUE 

XGRID=I SIZE/10 

I F ( MOD ( IS IZ E , 1 0 ) .NE.O) XGRI D=XGR I C + l 
CALL PLCTDV!X,Y, I S I Z E ,XGR I D , 3 , N 8 LK ) 
CONTINUE 



if ifififififififif 

PERFORM CORRECTIONS 

4: if if if if if if if # * if if if if if ^ if # # 



c 

c 

c 

c 

c 



DO 3010 1 = 1, N2 
FREQ= ( 1-1 )*SF/ISI ZE 
PI=3. 14159 

G=( !4.*PI*S IGMA *F REG /C)*C0S( THETA) )**2 
FINAL! I ) = FLOAT! I ARY! I ) ) 

I TEMP= I -I SF +1 

I TEMP2= I-I3M+1 , „ „„ 

IFCI.GE.ISF) FINAL! I ) =F I N AL ! I ) + EXP ( -G/ 2 . )*0*DS* 
SFINAL! I TEMP ) _ _ 

I F( I.GE .IBM ) FINAL! I )=FINAL ( I )-RCCEFF*D*DB* 
&FINAL! I TEMP2 ) 

3010 CONTINUE 

PRINT TIME PLOT AFTER CORRECTIONS 

££=(:* ******££;!;*## I*##:}: ####$* £* 4 =*#=!=* 



IF! IZ4.NE.YES) GOTO 3400 
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I F ( IZ5. NE.YES) GOTO 3350 
C 

C ***TIME PLOT FOR ENTIRE SIGNAL*** 

IF( INUM2.GT .1000) GOTO 3310 
NUM=I NU M2 
KPL0T S= 1 
GOTO 3320 

3310 KPL0TS=INUM2/1000+1 
NUM= 1000 

3320 DO 3340 I=1,KPL0TS 

I F ( ( I.EQ.KPL0TS) .AND. (I .NE.l) ) NUM = I NUM 2- ( K PLOT S- 1 ) * 
S1COO 

DC 3330 J = 1 , NU M 
Y ( J ) =F I NAL ( ( 1-1 )*1000+J) 

X < J ) = { ( I-1)*1000+J-1)/SF 
3330 CONTINUE 

XGRI0=NUM/1 0 

I F ( MODI NUM, 10) .NE .0) XGRID = XGRI C+l 
NBLK=- I 

CALL PLOTDV (X,Y ,NUM,XGRID,2 ,NBLK ) 

3340 CONTINUE 
GOTO 3400 

***TIME PLOT FOR INDIVIDUAL BLOCKS*** 

3350 DO 3370 1 = 1 ,KPTR 
NELK= I S TACK ( I ) 

DO 3360 J = 1 » I SI ZE 

Y(J) = FINAL( (NBLK-1 ) *1 SI ZE + J ) 

X( J)=( ( NBLK-1 J*I SI ZE+J-1 ) /SF 
3360 CONTINUE 

XGRID=ISIZF/10 

I F ( M0D( IS I Z E , 1 0 ) .NE.O) XGRI D=XGR I C+l 
CALL PLCTDVtX, Y, I S I Z E , XGR 1 0 , 2 , N E L K ) 

3370 CONTINUE 



PRINT FREQUENCY PLOTS 

34C0 I F ( IZ3. NE.YES) GOTO 3450 

***FREQUENC Y PLOT BEFORE CORRECTIONS*** 
DO 3430 1=1 , KPTR 
NBLK= ISTACK (I ) 

DC 3410 J = 1 , ISI ZE 

A ( J) = FLO AT ( I ARY ( ( NBLK-1 )* IS I ZE + J ) ) 

B ( J ) = 0 . 0 
3410 CONTINUE 

CALL CFFT2( A,B,ISIZE, I SIZE, ISI ZE,-1) 
ITEMP=I SI ZE/2 
DO 3420 J=1 , ITEMP 
Y ( J ) = A ( J ) **2+B ( J ) **2 
Y(J)=10*AL0G10(Y(J)) 

X( J)=( J-l )*SF/I SI ZE 
3420 CONTINUE 

XGRID=ITEMP/10 

IF (MODI ITEMP, 10) .NE.O) XGR ID=X GR IC + 1 
CALL PL OTDV(X,Y, ITEMP ,XGR I D , 1 , N B L K ) 

3430 CONTINUE 

***FREQUENC Y PLOT AFTER CORRECTIONS*** 
3450 I F( IZ2. NE.YES) GOTO 3500 
DO 3480 1=1, KPTR 
N8LK= 1ST ACK ( I ) 

DC 3460 J = 1 , 1 SI ZE 

A( J)=FI NAL ( (NBLK-1 )*ISIZE+J ) 

B ( J ) = 0 . 0 
346C CONTINUE 

CALL CFFT2( A,B, ISIZE, ISIZE, ISI ZE,-1) 

I TEMP=I SI ZE/2 
DO 3470 J=l, ITEMP 
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347C 

348C 

C 

C 

C 

C 

C 

3500 

3510 

C 

C 

c 

c 

c 

c 

2120 

2180 

2200 

2410 

2440 

2460 

2465 

2466 

2470 

2485 

2486 

2438 



Y ( J) =A( J)**2+B( J)**2 
Y (J)=10.*ALCG10(Y ( J ) ) 

X( J)=( J-l )*SF/I SI ZE 

CONTINUE 

XGRIO=I TEMP/10 

IF (MOD! ITEMPtlO) .NE.O) XGRI D=XGR I C + l 
CALL PLOTOV(X, Y, I TEMP ,XGR I D , 4 , N ELK ) 
CONTINUE 

■r t 'p'P t 'i* *r* 

CONCLUSION 
PRINT 3510 

FORMAT ( 1 ' t 'ARE YOU FINISHED? (Y/N)— •) 

READ 540, I Z 6 

IF( IZ6. EQ.YES) STOP 

IREAD=1 

GOTO 672 

DEBUG SUBCHK 

END 



^PLOT SUBROUTINE* 



SUBROUTINE PLOT DV ( X , Y , N , XGR I D, M , N 8 ) 

INTEGER D,XGRID,YGRID,AXIS 

DIMENSION Y(IOOO) , C7 ( 101 ) , 0 ( 6 ) , X ( 1000) , KA X I S ( 51 ) 

DATA IDASH/1H-/, I STAR/1H*/, ID0T/1H./ 

DATA IBAR/1HI/, IPLUS/1H+/, I BLANK/1 H /,IX/1HX/ 

AX I S=5 1 
YGR I D=6 
XGRID=XGR ID+1 
N 1=N-1 

Y 6= Y ( 1 ) 

Y1=Y( 1 ) 

DO 2200 1=1 , N 

IF(Y6-Y(I ) .GE.0.0) GOTO 2180 

Y 6 = Y ( I ) 

I F ( Yl-Y(I) .LE.0.0 ) GOTO 2200 
Y1=Y( I ) 

CGNT INUE 

S=Y1*(AXIS-1)/(Y6-Y1) 

X1=X( 1 ) 

XI 0=X( N ) 

0( 1)=Y1 

G ( 6 )=Y6 

I I X=XGR ID-1 

DO 2410 1 = 1, 1 1 X 

C7(I)=X((I-1)*1Q+1) 

CONTINUE 

C7(XGRID)=C7(XGRID-1 ) +10* ( X ( 2) - X ( 1 ) ) 

I F ( N. EQ . ( XGR ID-1 )*10) C7< XGR ID) =X( N) 

I I Y=YGR ID-1 
OG 2440 1=2, 1 1 Y 

0 ( I )= ( FLOAT ( I — 1 ) * (Y6-Y1 )/FLOAT (YGRID-1 ) )+Yl 

CONTINUE 

WRITE! 6,2460 

FORMAT (///,' •) 

IF(NB.GT.O) GOTO 2466 
Ne=-NB 

PRINT 2465, NB 

FORMAT ( *0' , 32X, • PLOT* , IX, 12 ) 

GOTO 2485 
PRINT 2470, NB 

FORMAT! 'O' ,32X, 'BLOCK’, IX, 12) 

I F ( M. E Q . 1 ) PRINT 2486 , _ _ „ „ ,_,. T T 

FORMAT! • ' , 17X , ' DB * ' S VS. FREQUENCY (BEFORE CORRECTION 

I FtM.EQ.2) PRINT 2488 

FORMAT ( * >, 20X, ' VOLTAGE VS. TIME (AFTER CORRECTIONS)’) 
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IF(M.EQ.3) PRINT 2487 

2487 FORMAT ( 1 1 , 20X , ' VOLT A GE VS. TIME (BEFORE COPR ECT I CNS ) ') 
IF(M.EQ.4) PRINT 2489 

2489 F C RM AT ( ' ',17X,'De"S VS. FREQUENCY (AFTER ', 

& • CORRECTIONS ) ' ) 

WRITE<6,2500) < 0 (I) , I =1 , YGR I D ) 

2500 F0RMAT(9X,11(1PE10.2) ) 

S 1= ( X10-X1 ) /10.0*(XGRID-1 ) 

D = 1 
L 1 = 1 
L = 1 

IZ=IFIX(-S+1.5) 

ITEMP= ( XGRID-1 )*10+1 
DO 2900 1 1 = 1 , ITE MP 
IFCN.LT.il) GOTO 2510 

YT EMP= ( Y ( II )*FLO AT (YGRID-1) *10.0/(Y6-Y1 ) )-S 
I Y=IF IX (YTEMP+1 • 5 ) 

2510 IF(Ll.GT.L) GOTO 2760 
DO 2650 I P= 1 , AXI S 
K AXI S ( I P ) = I DASH 
2650 CONTINUE 

DO 2680 1=1, AXIS, 10 
K AX I S ( I ) = I P LUS 
2680 CONTINUE 

IFCN.LT.il) GOTO 2720 

IF( (Yl.LE.O.O) .AND.(0.0.LE.Y6) ) KAXIS ( I Z ) = 1 COT 
K AX I S ( I Y ) = I STAR 

2720 WRITE(6 ,2725) C 7 ( D ) , ( KAX I S ( J ) , J = 1 , AX I S ) 

2725 FORMAT! 1PE12. 2, 2X,115A1) 

L1=L1+1 0 
D = D+1 
GCTO 2970 

2760 DC 2780 IP=1,AXIS 
K AXI S ( I P )= I BLANK 
2780 CONTINUE 

DC 2810 1 = 1, AXIS, 10 
KAXISt I ) = I B AR 
2810 CONTINUE 

IF(N.LT.Il) GOTO 2860 

IF({ Yl.LE.O.O). AN D.(0.0.LE.Y6)) KAXIS ( I Z ) = 1 COT 
K AX I S ( I Y ) = I STAR 

2860 WRITE (6, 2 86 5) ( KA X I S ( J ) , J = 1 , AX I S ) 

2865 FORMAT! 15X, 115A1 ) 

2870 L=L+1 

IF(L.GT .(XGRID-1 ) *10+2) GOTO 2910 
2900 CONTINUE 
2910 RETURN 

CEBUG SUBCHK 
END 
C 
C 

C ***FOUPIER TRANSFORM SUBROUTINE*** 

C 

SUBROUTINE CFFT2 ( A , B , NTOT , N , NSP AN , I SN ) 
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